Breeding Bird Surveys - Cleveland Metroparks (new and improved, with landscape covariates!)
Nathan W. Byer
2022-11-29
The dataset
Tim Krynak, Erik Shaffer, Jen Brumfield, and other Cleveland Metroparks and affiliated colleagues conducted Breeding Bird Surveys (BBS) at a total of 211 routes between 2017 and 2021. 164 of these routes overlapped with long-term vegetation monitoring efforts through the Cleveland Metroparks Plant Community Assessment Program (PCAP). This provided an excellent opportunity to explore the potential relationships between vegetative composition, landscape structure, and breeding bird ecology and community composition.
Initial processing
First, we just read in the full BBS dataset. For our purposes, we will remove any BBS observations not paired with PCAP plots. We can also visualize the distribution of survey effort through time.
CMP_bbs_pcap_join<-read.csv("G:/NaturalResources/Byer/datasets/BBS/CMP_bbs_pcapjoin_elimWC_EATO.csv",header=T)
CMP_bbs_pcap_join_pcaponly<-CMP_bbs_pcap_join[!is.na(CMP_bbs_pcap_join$Plot.ID),]
datetime<-strptime(CMP_bbs_pcap_join_pcaponly$survey_cre,format="%m/%d/%Y %H:%M")
hist(datetime,breaks="months")
pcap_1styear_landscape<-read.csv("G:/NaturalResources/Byer/datasets/BBS/pcap_4yr_landscape_variables_revLDI.csv",header=T)Next, we need to start grouping this dataset for downsteram analysis. This involves summarizing counts of each species for each plot. We also replace any NA values in the data frame with zeroes - this will be important for downstream applications, as we want to distinguish between unobserved species at a plot and check (typically coded 0s) and checks that were “skipped” (NAs).
CMP_bbs_pcap_join_pcaponly_pcapsum<-CMP_bbs_pcap_join_pcaponly %>%
group_by(simp_plot, reservatio,community,species_4_) %>%
summarise(n=n(),longitude=mean(xcoord), latitude=mean(ycoord),vibifq=mean(vibi_fq),vibifqnotrees=mean(vibi_fq_no_trees),carex=mean(carex_metric_value),cyper=mean(cyperaceae_metric_value),
dicot=mean(dicot_metric_value),shade=mean(shade_metric_value),shrub=mean(shrub_metric_value),
hydro=mean(hydrophyte_metric_value),svp=mean(svp_metric_value),ap=mean(ap_ratio_metric_value),
fqi=mean(fqai_score),bryo=mean(bryophyte_metric_value),hydro_rel=mean(hydrophyte_rel_cov_metric_value),
hydro_rel_notrees=mean(hydrophyte_rel_cov_no_trees),
sens_rel=mean(sensitive_rel_cov_metric_value),
sens_rel_notree=mean(sensitive_rel_cov_no_trees),
tol_rel=mean(tolerant_rel_cov_metric_value),
tol_rel_notree=mean(tolerant_rel_cov_no_trees),
inv=mean(invasive_graminoids_metric_value),
small_tree=mean(small_tree_metric_value),
subcanopy=mean(subcanopy_iv),
canopy=mean(canopy_iv.1)
) %>%
spread(species_4_,n) %>%
replace(is.na(.),0)
names(pcap_1styear_landscape)[names(pcap_1styear_landscape) == 'site.name'] <- 'simp_plot'
CMP_bbs_pcap_join_pcaponly_pcapsum_landscape<-left_join(x = CMP_bbs_pcap_join_pcaponly_pcapsum,y = pcap_1styear_landscape,by="simp_plot")
CMP_bbs_pcap_join_pcaponly_pcapsum<-column_to_rownames(CMP_bbs_pcap_join_pcaponly_pcapsum_landscape,var="simp_plot")This is a pretty unwieldly dataset at present - we have 164 rows (representing our BBS/PCAP plots) and 131 columns. Out of those 131 columns, 103 are unique species (designated with four-letter codes), so most of the first 27 covariates are essentially our “predictors” - in this case, our PCAP data.
Let’s do a bit more processing on this dataset to get it in shape for downstream analyses. First, we will separate out our bird community data.
birdcomm<-CMP_bbs_pcap_join_pcaponly_pcapsum[,28:(130)]
birdcomm[is.na(birdcomm)] <- 0Next, we separate out our PCAP data.
vegecomm<-CMP_bbs_pcap_join_pcaponly_pcapsum[,5:26]Finally, we will separate out the 1st year landscape covariates. As a quick note, these will not be strictly comparable with the 2nd survey landscape covariates, but should be pretty close (~5 years is likely not enough for substantial land cover change).
landcomm<-CMP_bbs_pcap_join_pcaponly_pcapsum[,c(142:152,156,157)]
landcomm[is.na(landcomm)] <- 0Since we are using the PCAP dataset as predictors here, we should likely check for correlations between these covariates. I have selected a (semi-)arbitrary cut-off of 0.6 to identify covariates as correlated; we then drop any covariates that are likely redundant.
vegecomm_cor<-cor(vegecomm)
set.seed(12345)
corvar<-findCorrelation(vegecomm_cor,cutoff=0.6)
vegecomm_drop<-vegecomm[,-corvar]We can also check for correlations between the landscape covariates.
cols<-seq(1,13,1)
landcomm[,cols] = apply(landcomm[,cols], 2, function(x) as.numeric(as.character(x)))
landcomm_cor<-cor(landcomm)
set.seed(12345)
corvar_land<-findCorrelation(landcomm_cor,cutoff=0.6)
landcomm_drop<-landcomm[,-corvar_land]We can then check correlations between PCAP and landscape covariates.
landvegcomm_drop<-cbind(vegecomm_drop,landcomm_drop)
landvegcomm_drop_cor<-cor(landvegcomm_drop)
set.seed(12345)
corvar_landveg<-findCorrelation(landvegcomm_drop_cor,cutoff=0.6)No covariates were overly correlated in this dataset, so we can proceed.
Now, we have:
- a vegecomm dataset - which representing PCAP measurements for each of the 164 BBS/PCAP plots.
- a birdcomm dataset - representing presence/absence of each bird species for each survey in each plot.
- a landcomm dataset - representing landscape covariates associated with these PCAP plots.
Note that we are not explicitly considering differences in bird counts (or occupancy, or detectability) through time at this point - we will get to that later. For now, we just want to explore broad associations between these three datasets!
Note as of 6/13/2022 - we will also want to only work with only binary data.
Community analyses
When working with datasets like these, we first have to acknowledge that these datasets are multivariate - with many possible predictor covariates (for example, plot-level characteristics) and many possible response covariates (for example, presence/absence of species in a plot). Fortunately, community ecology as a field has developed a number of standard approaches for exploring multivariate datasets, of which the most common is (arguably!) ordination. In essence, ordination techniques try to shrink down highly multidimensional datasets into a more manageable number of dimensions - often just a handful (and for initial plotting, typically just two). There are many, many possible ordination tools available for exploring community ecological datasets, and we will explore a few.
unconstrained ordination
In this section, we will be exploring ways of exploring the bird dataset on its own - basically, without any predictors.
First, we must reclassify our bird community data onto a binary scale - ranging from 0 (no detection) to 1 (detection).
We also set our color palette for residency status, so that Red = resident, Green = short-distance migrant, and Blue = neotropical migrant.
birdcomm_bin<-birdcomm
birdcomm_bin<- birdcomm_bin %>% mutate_if(is.numeric, ~1 * (. >= 1))
BBSMigrants<-read.csv("G:/NaturalResources/Byer/datasets/BBS/BBS Migants.csv",header=T)
BBSMigrants_sub<-BBSMigrants[BBSMigrants$SpeciesCode %in% colnames(birdcomm_bin),]
Residency<-BBSMigrants_sub$Residency
color_resid<-Residency
color_resid[color_resid %in% "NE"]<-"blue"
color_resid[color_resid %in% "M"]<-"green"
color_resid[color_resid %in% "R"]<-"red"
BBSMigrants_sub$color_resid<-color_residOne of the most common approaches for unconstrained ordination is called Principal Components Analysis (PCA). I have illustrated vegan’s rda() function (which we will use a bit more below!). Out of the available options, I find rda() to be much easier to use, so I tend to use it preferentially. Here is a plot with convex hulls and points colored by reservation.
while a bit hard to interpret, plots within Rocky River (RR) seem to be the most variable - and it might be that we see higher diversity in that area. We will see as we proceed if that is true.
metric multidimensional scaling (Principal Coordinates Analysis)
You can think of this as a distance-based analog to PCA. In this case, since we are using Bray-Curtis distances on binary data (presence/absence), this is equivalent to Sorensen’s index.
non-metric multidimensional scaling
bird species most associated with community differences
We can also look at which birds are most differentiated across communities for both the PCOA and NMDS ordinations. The below code essentially a) relates individual bird presence/absence to the final PCOA ordination, b) extracts p-values for associations with MDS1 and 2 (the first two ordination axes), and c) selects only those birds that are most differentiated at p less than or equal to 0.001. We will do so for the PCA, PCOA, and NMDS.
For all of these analyses, keep in mind that these are differentiated across plots, rather than reservations. We will explore a way to identify reservation-level differentiated species a litle bit below.
Here is the PCA:
## species PC1 PC2 pvals
## 1 ACFL -0.69716294 -0.1508935718 0.001
## 2 AMCR -0.24666198 -0.2963413159 0.001
## 3 AMGO 0.09382893 -0.2841531898 0.001
## 4 AMKE 0.38723531 -0.2455200453 0.001
## 5 BAOR 0.46434543 -0.1374142804 0.001
## 6 BARS 0.47425412 -0.2189233645 0.001
## 7 BCCH -0.25735224 -0.3604541492 0.001
## 8 BGGN -0.21605918 -0.2719538252 0.001
## 9 BHCO 0.17428135 -0.2790791107 0.001
## 10 BRTH 0.41973857 -0.0869533026 0.001
## 11 CANG 0.15587408 -0.2491690243 0.001
## 12 CEDW 0.50170158 -0.0806099475 0.001
## 13 COGR 0.45349980 0.1690696909 0.001
## 14 COYE 0.48298960 -0.2005910142 0.001
## 15 DEJU -0.35182512 -0.2234465042 0.001
## 16 EABL 0.21626153 -0.3348356997 0.001
## 17 EAKI 0.44209119 -0.2979537624 0.001
## 18 EAME 0.32949356 -0.4853999033 0.001
## 19 EAWP -0.52013724 -0.3050792447 0.001
## 20 EUST 0.64151056 -0.0626558506 0.001
## 21 FISP 0.25701439 -0.2476275133 0.001
## 22 GCFL -0.07114924 -0.3336146265 0.001
## 23 GRCA 0.51950507 -0.0464609980 0.001
## 24 HAWO -0.26479071 -0.3078359029 0.001
## 25 HOSP 0.32094904 0.1140673835 0.001
## 26 HOWA -0.44170756 -0.4240075369 0.001
## 27 HOWR 0.31383507 -0.1488100390 0.001
## 28 INBU 0.30394564 -0.3598407644 0.001
## 29 KILL 0.34044627 -0.2474971881 0.001
## 30 MODO 0.54350836 -0.2372128158 0.001
## 31 NOFL 0.32814765 -0.2185413306 0.001
## 32 NRWS 0.36829543 -0.0329679429 0.001
## 33 OROR 0.68797727 0.0475092454 0.001
## 34 PIWO -0.21798535 -0.2947734707 0.001
## 35 RBGR 0.29323586 -0.2267343583 0.001
## 36 REVI -0.38049218 -0.3368389029 0.001
## 37 RWBL 0.64080597 -0.0549021700 0.001
## 38 SCTA -0.53238990 -0.3546866667 0.001
## 39 SOSP 0.56532468 -0.1184872259 0.001
## 40 TRES 0.56126450 -0.1531979241 0.001
## 41 TUVU 0.24506841 -0.3094927503 0.001
## 42 WAVI 0.61917664 0.0134397815 0.001
## 43 WBNU -0.49381593 -0.2739331013 0.001
## 44 WIFL 0.47094574 0.0004034451 0.001
## 45 WOTH -0.31895854 -0.3625536846 0.001
## 46 YEWA 0.71017760 -0.0122296428 0.001
## 47 YTVI -0.28821198 -0.2674515592 0.00147 species were identified as significantly contributing to variation between plots. These are my best attempts to break these down by typical migratory patterns for each species in Ohio (courtesy of www.allaboutbirds.org).
Neotropical migrants: Acadian Fly-Catcher, Baltimore Oriole, Barn Swallow, Common Yellowthroat, Eastern Kingbird, Eastern Wood-Pewee, Great Crested Flycatcher, Hooded Warbler, Indigo Bunting, Orchard Oriole, Rose-Breasted Grosbeak, Red-Eyed Vireo, Scarlet Tanager, Warbling Vireo, Willow Flycatcher, Wood Thrush, Yellow Warbler, Yellow-Throated Vireo = 18 species
Resident: American Crow, American Kestrel, European Starling, Hairy Woodpecker, House Sparrow, Pileated Woodpecker, Song Sparrow, White-Breasted Nuthatch = 8 species
Short-distance migrant: American Goldfinch, Black-capped Chickadee, Blue-Gray Gnatcatcher, Brown-headed Cowbird, Brown Thrasher, Cedar Waxwing, Common Grackle, Eastern Bluebird, Eastern Meadowlark, Field Sparrow, Gray Catbird, House Wren, Killdeer, Mourning Dove, Northern Flicker, Northern Rough-winged Swallow, Red-winged Blackbird, Tree Swallow, Turkey Vulture = 19 species
Resident or migratory: Canada Goose, Dark-eyed Junco = 2 species
We can repeat this procedure for the PCOA results too. Whereas PCA (above) was based on summarizing our bird community in the fewest number of axes that explain the most variance, PCOA is trying to maximize distance between plots. This is somewhat semantic in many cases, but will lead to differing representations of which birds are most associated with ordination axes.
Anyway! on to the PCOA:
## species MDS1 MDS2 pvals
## 1 ACFL -0.731621770 -0.21657507 0.001
## 2 AMRO 0.212185897 0.26237816 0.001
## 3 BAOR 0.496278108 -0.11707782 0.001
## 4 BARS 0.396414060 -0.14525664 0.001
## 5 BGGN -0.253771660 -0.21321227 0.001
## 6 BHCO 0.252189128 -0.33826670 0.001
## 7 BRTH 0.332695547 -0.17755430 0.001
## 8 CEDW 0.528291870 -0.15714239 0.001
## 9 CHSW -0.008700066 -0.32221685 0.001
## 10 COGR 0.449821060 0.15723190 0.001
## 11 COYE 0.516757753 -0.36308017 0.001
## 12 DEJU -0.353475491 -0.17838667 0.001
## 13 EAKI 0.323351441 -0.09112853 0.001
## 14 EATO 0.130826116 -0.28201583 0.001
## 15 EAWP -0.531425079 -0.06810441 0.001
## 16 EUST 0.564199285 0.00876616 0.001
## 17 FISP 0.247783137 -0.21620323 0.001
## 18 GRCA 0.606185199 -0.08760151 0.001
## 19 HAWO -0.240271017 -0.22737147 0.001
## 20 HOFI 0.344663351 0.14118008 0.001
## 21 HOSP 0.357013239 0.15621096 0.001
## 22 HOWA -0.446954750 -0.49945770 0.001
## 23 HOWR 0.363416944 0.14393149 0.001
## 24 INBU 0.358710671 -0.48692796 0.001
## 25 LOWA -0.216453442 -0.34516208 0.001
## 26 MODO 0.503095502 -0.13966980 0.001
## 27 NOFL 0.338135730 -0.04580320 0.001
## 28 NRWS 0.299804802 -0.05498756 0.001
## 29 OROR 0.635092702 -0.04190854 0.001
## 30 OVEN -0.178146781 -0.19698291 0.001
## 31 PIWO -0.211414515 -0.44190745 0.001
## 32 RBGR 0.335449325 -0.27805077 0.001
## 33 REVI -0.427474339 -0.05322856 0.001
## 34 RWBL 0.659017118 -0.04422247 0.001
## 35 SCTA -0.563727972 -0.35767304 0.001
## 36 SOSP 0.654370712 -0.10449673 0.001
## 37 TRES 0.446180461 -0.04558460 0.001
## 38 VEER -0.141621359 -0.32912045 0.001
## 39 WAVI 0.590289204 0.02894698 0.001
## 40 WBNU -0.423708109 -0.15390376 0.001
## 41 WIFL 0.417164492 -0.19105681 0.001
## 42 WOTH -0.273041942 -0.51805873 0.001
## 43 YEWA 0.695495064 -0.12652046 0.001
## 44 YTVI -0.257925903 -0.31119268 0.001
So, now we have 44 species that contribute most to community distance between plots.
Neotropical migrants: Acadian Fly-Catcher, Baltimore Oriole, Barn Swallow, Chimney Swift, Common Yellowthroat, Eastern Kingbird, Eastern Wood-Pewee, Hooded Warbler, Indigo Bunting, Louisiana Water-Thrush, Orchard Oriole, Ovenbird, Rose-Breasted Grosbeak, Red-Eyed Vireo, Scarlet Tanager, Veery, Warbling Vireo, Willow Flycatcher, Wood Thrush, Yellow Warbler, Yellow-Throated Vireo = 21 species
Resident: European Starling, Hairy Woodpecker, House Finch, House Sparrow, Pileated Woodpecker, Song Sparrow, White-Breasted Nuthatch = 7 species
Short-distance migrant: American Robin, Blue-Gray Gnatcatcher, Brown-headed Cowbird, Brown Thrasher, Cedar Waxwing, Common Grackle, Eastern Towhee, Field Sparrow, Gray Catbird, House Wren, Mourning Dove, Northern Flicker, Northern Rough-winged Swallow, Red-winged Blackbird, Tree Swallow = 15 species
Resident or migratory: Dark-eyed Junco = 1 species
and now the NMDS:
## species NMDS1 NMDS2 pvals
## 1 ACFL -0.6618348 -0.341392498 0.001
## 2 AMCR -0.1809373 -0.228515977 0.001
## 3 AMRO 0.1629920 0.302705875 0.001
## 4 BAOR 0.4772962 -0.127988404 0.001
## 5 BARS 0.3860442 -0.157950082 0.001
## 6 BCCH -0.1840405 -0.247342387 0.001
## 7 BGGN -0.1661447 -0.310696158 0.001
## 8 BHCO 0.2053812 -0.323500671 0.001
## 9 BRTH 0.2829251 -0.255309288 0.001
## 10 CEDW 0.4945065 -0.161253661 0.001
## 11 COGR 0.4097261 0.223524723 0.001
## 12 COYE 0.4294219 -0.454069222 0.001
## 13 DEJU -0.2800585 -0.283832953 0.001
## 14 EAKI 0.3026501 -0.153678817 0.001
## 15 EAWP -0.5238803 -0.171371275 0.001
## 16 EUST 0.5579432 0.034119785 0.001
## 17 FISP 0.1842373 -0.297461046 0.001
## 18 GRCA 0.5913728 0.010163788 0.001
## 19 HOFI 0.2689176 0.261298811 0.001
## 20 HOSP 0.2863387 0.281273735 0.001
## 21 HOWA -0.3076686 -0.569869867 0.001
## 22 HOWR 0.3001183 0.252110188 0.001
## 23 INBU 0.2709295 -0.507591667 0.001
## 24 LOWA -0.1380231 -0.364701465 0.001
## 25 MODO 0.4666425 -0.213629459 0.001
## 26 NOFL 0.3201925 0.012743833 0.001
## 27 OROR 0.6399450 -0.052994011 0.001
## 28 PIWO -0.1194412 -0.400313818 0.001
## 29 RBGR 0.2673543 -0.328081328 0.001
## 30 REVI -0.4342023 -0.084154867 0.001
## 31 RWBL 0.6464910 -0.020637450 0.001
## 32 SCTA -0.4227056 -0.474345456 0.001
## 33 SOSP 0.6408078 -0.111264386 0.001
## 34 TRES 0.4323321 -0.091271373 0.001
## 35 WAVI 0.5938069 -0.002640012 0.001
## 36 WBNU -0.3459764 -0.293921549 0.001
## 37 WIFL 0.3889553 -0.241823500 0.001
## 38 WOTH -0.1553185 -0.520658575 0.001
## 39 YEWA 0.6794610 -0.197506714 0.001
## 40 YTVI -0.1666356 -0.363636777 0.001So, now we have 40 species that contribute most to community distance between plots.
Neotropical migrants: Acadian Fly-Catcher, Baltimore Oriole, Barn Swallow, Common Yellowthroat, Eastern Kingbird, Eastern Wood-Pewee, Hooded Warbler, Indigo Bunting, Louisiana Water-Thrush, Orchard Oriole, Rose-Breasted Grosbeak, Red-Eyed Vireo, Scarlet Tanager, Warbling Vireo, Willow Flycatcher, Wood Thrush, Yellow Warbler, Yellow-Throated Vireo = 18 species
Resident: American Crow, European Starling, House Finch, House Sparrow, Pileated Woodpecker, Song Sparrow, White-Breasted Nuthatch = 7 species
Short-distance migrant: American Robin, Black-Capped Chickadee, Blue-Gray Gnatcatcher, Brown-headed Cowbird, Brown Thrasher, Cedar Waxwing, Common Grackle, Field Sparrow, Gray Catbird, House Wren, Mourning Dove, Northern Flicker, Red-winged Blackbird, Tree Swallow = 14 species
Resident or migratory: Dark-eyed Junco = 1 species
While these three sets of species are pretty similar, we can also take a look at which species are different across approaches.
outersect(fit_spp$species,fit_pcoa_spp$species,fitpca_spp$species)## [1] "CHSW" "EATO" "OVEN" "VEER" "AMGO" "AMKE" "CANG" "EABL" "EAME" "GCFL" "KILL" "TUVU"So American Goldfinch, American Kestrel, Canada Goose, Chimney Swift, Eastern Bluebird, Eastern Meadowlark, Easter Towhee, Great-Crested Flycatcher, Killdeer, Ovenbird, Tufted Titmouse, Turkey Vulture, and Veery. At least to me, nothing particularly sticks out as a commonality across these species, but it may make just as much sense to simply pool these species across approaches:
sort(unique(c(fit_spp$species,fit_pcoa_spp$species,fitpca_spp$species)))## [1] "ACFL" "AMCR" "AMGO" "AMKE" "AMRO" "BAOR" "BARS" "BCCH" "BGGN" "BHCO" "BRTH" "CANG" "CEDW" "CHSW" "COGR" "COYE" "DEJU" "EABL" "EAKI" "EAME" "EATO" "EAWP" "EUST" "FISP" "GCFL" "GRCA" "HAWO" "HOFI" "HOSP" "HOWA" "HOWR" "INBU" "KILL" "LOWA" "MODO" "NOFL" "NRWS" "OROR" "OVEN" "PIWO" "RBGR" "REVI" "RWBL" "SCTA" "SOSP" "TRES" "TUVU" "VEER" "WAVI" "WBNU" "WIFL" "WOTH" "YEWA" "YTVI"This gives us a list of 54 bird species that seem to be particularly differentiated/variable across plots. That is about half of the total number of species in the overall dataset.
As we proceed, we will begin layering on our environmental datasets using three main approaches: constrained ordination, linear models of species richness, and multi-species occupancy models.
constrained
Whereas the ordination approaches above were focused on unconstrained ordination - in this case, exploring how species are arranged between sample sites - there is an entire suite of approaches that address how environmental gradients shape community assemblages among sample sites. These are called constrained ordination approaches. I am personally a huge proponent of redundancy analysis (RDA) as a flexible constrained ordination approach. Taken from the GUSTA ME tutorial on RDA:
“RDA can also be considered a constrained version of principal components analysis (PCA), wherein canonical axes - built from linear combinations of response variables - must also be linear combinations of the explanatory variables (i.e. fitted by [multiple linear regression]). The RDA approach generates one ordination in the space defined by the matrix of response variables and another in the space defined by the matrix of explanatory variables.”
So, in essence, the ordination plot produced will depict an optimal arrangement of environmental axes (RDA1, 2, etc.) that maximally explains variation in the species dataset.
We are going to use stepwise selection to explore which PCAP covariates most explain plot-level differences.
p<-rda(birdcomm_bin~.,vegecomm_drop,scale=T)
fit<-envfit(p,vegecomm_drop,scale=T)
p_0<-rda(birdcomm_bin~1,vegecomm_drop,scale=T)
p_1<-rda(birdcomm_bin~.,vegecomm_drop,scale=T)
set.seed(12345)
p_step<- ordistep(p_0, scope = formula(p_1),trace=F)
p_step## Call: rda(formula = birdcomm_bin ~ tol_rel + dicot + cyper + svp, data = vegecomm_drop, scale = T)
##
## Inertia Proportion Rank
## Total 103.00000 1.00000
## Constrained 7.68058 0.07457 4
## Unconstrained 95.31942 0.92543 102
## Inertia is correlations
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3 RDA4
## 4.484 1.491 0.984 0.721
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 6.460 3.678 3.582 3.215 3.037 2.621 2.553 2.371
## (Showing 8 of 102 unconstrained eigenvalues)anova(p_step)## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = birdcomm_bin ~ tol_rel + dicot + cyper + svp, data = vegecomm_drop, scale = T)
## Df Variance F Pr(>F)
## Model 4 7.681 3.2029 0.001 ***
## Residual 159 95.319
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(p_step,by="term")## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = birdcomm_bin ~ tol_rel + dicot + cyper + svp, data = vegecomm_drop, scale = T)
## Df Variance F Pr(>F)
## tol_rel 1 4.048 6.7529 0.001 ***
## dicot 1 1.362 2.2727 0.001 ***
## cyper 1 1.179 1.9674 0.002 **
## svp 1 1.090 1.8188 0.004 **
## Residual 159 95.319
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(p_step,by="axis")## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = birdcomm_bin ~ tol_rel + dicot + cyper + svp, data = vegecomm_drop, scale = T)
## Df Variance F Pr(>F)
## RDA1 1 4.484 7.4791 0.001 ***
## RDA2 1 1.491 2.4876 0.001 ***
## RDA3 1 0.984 1.6421 0.007 **
## RDA4 1 0.721 1.2030 0.120
## Residual 159 95.319
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1pcapplot<-ordiplot(p_step,scaling=3,type="n")
points(pcapplot, "sites",pch=3, col=colpal_reserv[reserv])
ordihull(pcapplot,groups = reserv, col=colpal_reserv, lwd=3)
title("PCA (vegan rda version)")
legend(x="right", legend=levels(reserv), col=colpal_reserv,pch = 3,cex=0.5)
ggord(p_step, xlim=c(-1.25,1),ylim = c(-0.5, 1),txt=4,arrow=TRUE,ptslab=TRUE,addsize=3, size =4,addcol = color_resid)
So this tells us that the most explanatory group of PCAP covariates is relative cover of tolerant plants (tol_rel), cover of dicots, cover of seedless vascular plants (svp), and cover of sedges (cyper). Rather conveniently, the relative cover of dicots, seedless vascular plants, and sedges is negatively associated with the relative cover of tolerant plant species.
As a quick caveat, these four variables only explain about 7.3% of the variation, as seen here (red is our environmental axes, black is our unconstrained axes):
constrained_eig <- p_step$CCA$eig/p_step$tot.chi*100
unconstrained_eig <- p_step$CA$eig/p_step$tot.chi*100
expl_var <- c(constrained_eig, unconstrained_eig)
barplot (expl_var[1:20], col = c(rep ('red', length (constrained_eig)), rep ('black', length (unconstrained_eig))),
las = 2, ylab = '% variation')
This is not that terribly uncommon for redundancy analyses - It is often difficult to find environmental gradients that explain all of the variance in large datasets like ours. Biologically, we likely do not necessarily expect that a single environmental gradient alone drives all of the variation in bird presence/absence, and intuitively know that we also have to worry about observer effects, weather, and other unmeasured characteristics. Regardless of any of the above, the ordination is significant - as are the relationships with these covariates, so we can safely say that nearly 10% of the variation in our bird species dataset is explainable by these 4 PCAP covariates!
We can also explore just the landscape covariate associations:
p_land<-rda(birdcomm_bin~.,landcomm_drop,scale=T)
fit<-envfit(p_land,landcomm_drop,scale=T)
p_land_0<-rda(birdcomm_bin~1,landcomm_drop,scale=T)
p_land_1<-rda(birdcomm_bin~.,landcomm_drop,scale=T)
set.seed(12345)
p_land_step<- ordistep(p_land_0, scope = formula(p_land_1),trace=F)
anova(p_land_step)## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = birdcomm_bin ~ ldi_3k + dist_road + activity.area + developed.edge, data = landcomm_drop, scale = T)
## Df Variance F Pr(>F)
## Model 4 5.901 2.4157 0.001 ***
## Residual 159 97.099
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(p_land_step,by="term")## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = birdcomm_bin ~ ldi_3k + dist_road + activity.area + developed.edge, data = landcomm_drop, scale = T)
## Df Variance F Pr(>F)
## ldi_3k 1 2.572 4.2124 0.001 ***
## dist_road 1 1.483 2.4281 0.001 ***
## activity.area 1 0.967 1.5840 0.019 *
## developed.edge 1 0.878 1.4384 0.033 *
## Residual 159 97.099
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(p_land_step,by="axis")## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = birdcomm_bin ~ ldi_3k + dist_road + activity.area + developed.edge, data = landcomm_drop, scale = T)
## Df Variance F Pr(>F)
## RDA1 1 3.655 5.9847 0.001 ***
## RDA2 1 0.880 1.4418 0.207
## RDA3 1 0.773 1.2659 0.265
## RDA4 1 0.593 0.9707 0.507
## Residual 159 97.099
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ordiplot(p_land_step,scaling=3)
constrained_eig <- p_land_step$CCA$eig/p_land_step$tot.chi*100
unconstrained_eig <- p_land_step$CA$eig/p_land_step$tot.chi*100
expl_var <- c(constrained_eig, unconstrained_eig)
barplot (expl_var[1:20], col = c(rep ('red', length (constrained_eig)), rep ('black', length (unconstrained_eig))),
las = 2, ylab = '% variation')
And now for the vegetation and land cover dataset:
p_landveg<-rda(birdcomm_bin~.,landvegcomm_drop,scale=T)
fit_landveg<-envfit(p_landveg,landvegcomm_drop,scale=T)
p_landveg_0<-rda(birdcomm_bin~1,landvegcomm_drop,scale=T)
p_landveg_1<-rda(birdcomm_bin~.,landvegcomm_drop,scale=T)
set.seed(12345)
p_landveg_step<- ordistep(p_landveg_0, scope = formula(p_landveg_1),trace = F)
p_landveg_step## Call: rda(formula = birdcomm_bin ~ tol_rel + ldi_3k + dist_road + dicot + cyper + svp + developed.edge, data = landvegcomm_drop, scale = T)
##
## Inertia Proportion Rank
## Total 103.0000 1.0000
## Constrained 11.4155 0.1108 7
## Unconstrained 91.5845 0.8892 102
## Inertia is correlations
##
## Eigenvalues for constrained axes:
## RDA1 RDA2 RDA3 RDA4 RDA5 RDA6 RDA7
## 5.505 1.791 1.420 0.877 0.755 0.579 0.489
##
## Eigenvalues for unconstrained axes:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
## 5.867 3.544 3.431 3.156 2.597 2.507 2.357 2.328
## (Showing 8 of 102 unconstrained eigenvalues)anova(p_landveg_step)## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = birdcomm_bin ~ tol_rel + ldi_3k + dist_road + dicot + cyper + svp + developed.edge, data = landvegcomm_drop, scale = T)
## Df Variance F Pr(>F)
## Model 7 11.415 2.7778 0.001 ***
## Residual 156 91.585
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(p_landveg_step,by="term")## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = birdcomm_bin ~ tol_rel + ldi_3k + dist_road + dicot + cyper + svp + developed.edge, data = landvegcomm_drop, scale = T)
## Df Variance F Pr(>F)
## tol_rel 1 4.048 6.8957 0.001 ***
## ldi_3k 1 2.314 3.9410 0.001 ***
## dist_road 1 1.284 2.1872 0.002 **
## dicot 1 1.072 1.8267 0.003 **
## cyper 1 1.048 1.7854 0.002 **
## svp 1 0.849 1.4464 0.024 *
## developed.edge 1 0.800 1.3620 0.064 .
## Residual 156 91.585
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(p_landveg_step,by="axis")## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
##
## Model: rda(formula = birdcomm_bin ~ tol_rel + ldi_3k + dist_road + dicot + cyper + svp + developed.edge, data = landvegcomm_drop, scale = T)
## Df Variance F Pr(>F)
## RDA1 1 5.505 9.3762 0.001 ***
## RDA2 1 1.791 3.0501 0.001 ***
## RDA3 1 1.420 2.4196 0.001 ***
## RDA4 1 0.877 1.4934 0.180
## RDA5 1 0.755 1.2861 0.380
## RDA6 1 0.579 0.9867 0.844
## RDA7 1 0.489 0.8324 0.840
## Residual 156 91.585
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ordiplot(p_landveg_step,scaling=3)
constrained_eig <- p_landveg_step$CCA$eig/p_landveg_step$tot.chi*100
unconstrained_eig <- p_landveg_step$CA$eig/p_landveg_step$tot.chi*100
expl_var <- c(constrained_eig, unconstrained_eig)
barplot (expl_var[1:20], col = c(rep ('red', length (constrained_eig)), rep ('black', length (unconstrained_eig))),
las = 2, ylab = '% variation')
Some interesting things are going on here. RDA1 primarily captures both plot-scale and landscape-scale indicators of disturbance - namely, relative cover of tolerant plant species (plot level) and landscape development index (at a landscape-level) contrasted with indicators of less disturbance - namely, further distances to developed land, more seedless vascular plants, more cyperaceae, and more dicots. This one axis explains upwards of 5% of variation - so we are doing pretty well in that regard. Across all seven axes, we are approaching about 11% explained variance. In addition, stepwise selection basically confirms variable selection for the land cover and vegetation datasets run separately - the same variables included in those best-supported model statements are in this full model.
We can also explore bird species most significantly associated with our constrained ordination axes, too. We will focus on the full (land cover and vegetation) ordination.
# envfit() takes the output of metaMDS() and the species matrix you created
fitconcrda <- envfit(p_landveg_step, birdcomm_bin, perm = 999)
# extract p-values for each species
fitconcrda_pvals <- fitconcrda$vectors$pvals %>%
as.data.frame() %>%
rownames_to_column("species") %>%
dplyr::rename("pvals" = ".")
# extract coordinates for species, only keep species with p-val = 0.001
fitconcrda_spp <- fitconcrda %>%
scores(., display = "vectors") %>%
as.data.frame() %>%
rownames_to_column("species") %>%
full_join(., fitconcrda_pvals, by = "species") %>%
filter(pvals == 0.001)OK - this is getting towards our main analytical goals! Let’s sort by scores for that first environmental axis and plot these in order:
BBSMigrants_sub_sig<-BBSMigrants_sub[BBSMigrants_sub$SpeciesCode %in% fitconcrda_spp$species,]
BBSMigrants_sub_sig_reorder<-BBSMigrants_sub_sig[order(BBSMigrants_sub_sig$SpeciesCode),]
fitconcrda_spp_reorder<-fitconcrda_spp[order(fitconcrda_spp$species),]
fitconcrda_spp_join<-cbind(fitconcrda_spp_reorder,BBSMigrants_sub_sig_reorder)
fitconcrda_spp_order<-fitconcrda_spp_join[order(fitconcrda_spp_join$RDA1),]
fitconcrda_spp_order$species<-factor(fitconcrda_spp_order$species,levels=fitconcrda_spp_order$species)
fitconcrda_spp_order$color_resid[fitconcrda_spp_order$Residency %in% "NE"]<-"blue"
fitconcrda_spp_order$color_resid[fitconcrda_spp_order$Residency %in% "M"]<-"green"
fitconcrda_spp_order$color_resid[fitconcrda_spp_order$Residency %in% "R"]<-"red"
color_sub<-fitconcrda_spp_order$color_resid
ggplot(aes(x=RDA1,y=species),data=fitconcrda_spp_order)+geom_point()+theme_classic(base_size=18)+theme(axis.text.y = element_text(colour = color_sub))## Warning: Vectorized input to `element_text()` is not officially supported.
## ℹ Results may be unexpected or may change in future versions of ggplot2.
This is a pretty neat plot - Acadian Flycatcher, Scarlet Tanager, Eastern Wood-Pewee, White-breasted Nuthatch, Hooded Warbler, Red-Eyed Vireo, Dark-Eyed Junco, Wood Thrush, Hairy Woodpeckers, and Lousiana Water-Thrush all appear to be associated with higher sedge, seedless vascular plant, and dicot cover and lower cover of tolerant plants. In contrast, Yellow Warblers, Orchard Orioles, Red-Winged Blackbirds, European Starlings, Warbling Vireos, Song Sparrows, Grey Catbirds, Tree Swallows, Mourning Doves, and Cedar Waxwings are all associated with more tolerant plants. In addition, virtually all of our more tolerant vegetation-associated species are short-distance migrants, whereas the less tolerant vegetation-associated species are mostly neotropical migrants.
constrained - sites
While we have already explored environmental associations, we might also want to know what species are most differentiated between reservations. We will then use similar code to that used above to explore species that are most differentiated between reservations.
p_reservations<-rda(birdcomm_bin~reservatio,CMP_bbs_pcap_join_pcaponly_pcapsum,scale=T)
reservations_fitbirds<-envfit(p_reservations,birdcomm_bin)
reservations_fit<-envfit(p_reservations~reservatio,data=CMP_bbs_pcap_join_pcaponly_pcapsum,scale=T)
CMP_bbs_pcap_join_pcaponly_pcapsum$site<-rownames(CMP_bbs_pcap_join_pcaponly_pcapsum)
plot_df_reservations <- scores(p_reservations, display = "sites") %>%
as.data.frame() %>%
rownames_to_column("site") %>%
full_join(CMP_bbs_pcap_join_pcaponly_pcapsum, by = "site")
plot_rda_reservations <- ggplot(plot_df_reservations, aes(x = RDA1, y = RDA2, color = reservatio)) +
geom_point(size = 3, alpha = 0.8) +
stat_ellipse(linetype = 2, size = 1) +
labs(title = "RDA")
plot_rda_reservations
fitreservations_pvals <- reservations_fitbirds$vectors$pvals %>%
as.data.frame() %>%
rownames_to_column("species") %>%
dplyr::rename("pvals" = ".")
# extract coordinates for species, only keep species with p-val = 0.001
fit_reservations_spp <- reservations_fitbirds %>%
scores(., display = "vectors") %>%
as.data.frame() %>%
rownames_to_column("species") %>%
full_join(., fitreservations_pvals, by = "species") %>%
filter(pvals == 0.001)
fit_reservations_spp## species RDA1 RDA2 pvals
## 1 ACFL -0.66777571 0.117759129 0.001
## 2 AMCR -0.38752464 -0.218120450 0.001
## 3 AMKE 0.36069599 0.005207019 0.001
## 4 AMRE -0.11691782 -0.417572542 0.001
## 5 AMRO 0.26690627 0.242961322 0.001
## 6 BAOR 0.40912220 0.147084861 0.001
## 7 BARS 0.36426344 -0.297948461 0.001
## 8 BCCH -0.25654264 0.173591143 0.001
## 9 BRCR -0.28687748 -0.100539363 0.001
## 10 BRTH 0.24661825 -0.249101923 0.001
## 11 CEDW 0.40587243 -0.318529350 0.001
## 12 COGR 0.51500011 0.167583353 0.001
## 13 COYE 0.09667555 -0.532851862 0.001
## 14 DEJU -0.37639933 -0.016388874 0.001
## 15 EATO -0.13543382 -0.284272455 0.001
## 16 EAWP -0.46229095 0.436908082 0.001
## 17 EUST 0.54673983 -0.420762254 0.001
## 18 FISP -0.03106708 -0.366765924 0.001
## 19 GRCA 0.46354706 -0.107291445 0.001
## 20 HOFI 0.37928897 -0.073545643 0.001
## 21 HOSP 0.47178576 0.090612801 0.001
## 22 HOWA -0.57041850 -0.263142828 0.001
## 23 HOWR 0.31605596 0.122088141 0.001
## 24 INBU 0.14221796 -0.277445296 0.001
## 25 LOWA -0.30349510 -0.265588533 0.001
## 26 MALL 0.36539436 0.040174144 0.001
## 27 MODO 0.35803008 -0.169949443 0.001
## 28 NOFL 0.30190291 0.142456672 0.001
## 29 NRWS 0.44565648 0.077961759 0.001
## 30 OROR 0.53857451 -0.399968564 0.001
## 31 OVEN -0.32654527 -0.281863336 0.001
## 32 RBGU 0.36628682 -0.086059767 0.001
## 33 REVI -0.30390875 0.295754816 0.001
## 34 RWBL 0.50404846 -0.271167918 0.001
## 35 SCTA -0.54375433 0.062745434 0.001
## 36 SOSP 0.45897351 -0.181014710 0.001
## 37 TRES 0.44363253 -0.106280777 0.001
## 38 TUTI -0.28142388 0.201253917 0.001
## 39 TUVU 0.05394872 -0.299642895 0.001
## 40 VEER -0.25657805 -0.198795371 0.001
## 41 WAVI 0.54661431 -0.284150723 0.001
## 42 WBNU -0.47299945 0.272560858 0.001
## 43 WIFL 0.30162601 -0.294272609 0.001
## 44 WOTH -0.36226821 -0.032288378 0.001
## 45 YEWA 0.50814724 -0.424006710 0.001
## 46 YTVI -0.36532038 -0.235368151 0.001fit_reservations <- reservations_fit %>%
scores(., display = "factors") %>%
as.data.frame() %>%
rownames_to_column("reservations")
rda_reservation_plot_new <- ggplot(plot_df_reservations, aes(x = RDA1, y = RDA2)) +coord_fixed() +
geom_point(aes(color = reservatio), size = 3, alpha = 0.8) +
geom_segment(data = fit_reservations_spp, aes(x = 0, xend = RDA1, y = 0, yend = RDA2),col = "black") +
geom_text(data = fit_reservations_spp, aes(label = species))+
geom_text(data=as.data.frame(fit_reservations),aes(label=reservations))
rda_reservation_plot_new
Interesting! So when we consider reservation-of-origin, Hinckley seems to have a pretty unique bird community overall, and groups negatively along RDA1. In contrast, Brookside, Ohio and Erie Canal, and Washington all group positively along RDA1 - and are all quite heavily urbanized.
While I personally think the environmental gradient analyses are a bit more interesting, this does highlight reservations that are more or less distinct in bird community space. We can also plot just the reservations alone too:
ggplot(plot_df_reservations, aes(x = RDA1, y = RDA2)) +
coord_fixed() +
geom_text(data=as.data.frame(fit_reservations),aes(label=reservations))+xlim(-2,4)+ylim(-4,2)
lda of reservations
reservation<-factor(CMP_bbs_pcap_join_pcaponly_pcapsum$reservatio)
reservation_min3<-reservation[!(reservation %in% c("BK","EC","GP","HU","WA"))]
birdcomm_bin_min3<-birdcomm_bin[!(CMP_bbs_pcap_join_pcaponly_pcapsum$reservatio %in% c("BK","EC","GP","HU","WA")),]
birdcomm_bin_min3_drop1<-birdcomm_bin_min3[colSums(birdcomm_bin_min3) > 1]
color_resid_dropconstant<-color_resid[(colnames(birdcomm_bin) %in% colnames(birdcomm_bin_min3_drop1))]
lda_birdreserv<-lda(birdcomm_bin_min3_drop1, grouping = reservation_min3)
ggord(lda_birdreserv, reservation_min3, xlim=c(-15,15),ylim = c(-20, 20),txt=NULL,arrow=NULL)
ggord(lda_birdreserv, reservation_min3, xlim=c(-25,25),ylim = c(-20, 20),txt=NULL,arrow=NULL,alpha=0,alpha_el=0.5)
ggord(lda_birdreserv, reservation_min3, xlim=c(-25,25),ylim = c(-20, 20),alpha=0,alpha_el=0.5,veccol=color_resid_dropconstant,labcol=color_resid_dropconstant)
ggord(lda_birdreserv, reservation_min3, xlim=c(-25,25),ylim = c(-20, 20),veccol=color_resid_dropconstant,labcol=color_resid_dropconstant)
ggord(lda_birdreserv, reservation_min3, xlim=c(-10,10),ylim = c(-10, 10),alpha=0,alpha_el=0.5,veccol=color_resid_dropconstant,labcol=color_resid_dropconstant)
ggord(lda_birdreserv, reservation_min3, xlim=c(-7.5,7.5),ylim = c(-5, 6),alpha=0,alpha_el=0.5,var_sub = c("ACFL","HOWA","PIWO","SOSP"))
spe.class <- predict(lda_birdreserv)$class
# Posterior probabilities that the objects belong to those
# groups
spe.post <- predict(lda_birdreserv)$posterior
# Table of prior vs. predicted classifications
(spe.table <- table(reservation_min3, spe.class))
# Proportion of corrected classification
diag(prop.table(spe.table, 1))Community associations
Now, we can also look at how bird communities are differentiated between community types.
community<-factor(CMP_bbs_pcap_join_pcaponly_pcapsum$community.x)
p_community<-rda(birdcomm_bin~community,CMP_bbs_pcap_join_pcaponly_pcapsum,scale=T)
community_fitbirds<-envfit(p_community,birdcomm_bin)
community_fit<-envfit(p_community~community,data=CMP_bbs_pcap_join_pcaponly_pcapsum,scale=T)
CMP_bbs_pcap_join_pcaponly_pcapsum$site<-rownames(CMP_bbs_pcap_join_pcaponly_pcapsum)
plot_df_community <- scores(p_community, display = "sites") %>%
as.data.frame() %>%
rownames_to_column("site") %>%
full_join(CMP_bbs_pcap_join_pcaponly_pcapsum, by = "site")
plot_rda_community <- ggplot(plot_df_community, aes(x = RDA1, y = RDA2, color = community)) +
geom_point(size = 3, alpha = 0.8) +
stat_ellipse(linetype = 2, size = 1) +
labs(title = "RDA")
plot_rda_community
fitcommunity_pvals <- community_fitbirds$vectors$pvals %>%
as.data.frame() %>%
rownames_to_column("species") %>%
dplyr::rename("pvals" = ".")
# extract coordinates for species, only keep species with p-val = 0.001
fit_community_spp <- community_fitbirds %>%
scores(., display = "vectors") %>%
as.data.frame() %>%
rownames_to_column("species") %>%
full_join(., fitcommunity_pvals, by = "species") %>%
filter(pvals == 0.001)
fit_community_spp## species RDA1 RDA2 pvals
## 1 ACFL -0.70195248 0.097941485 0.001
## 2 AMCR -0.26483267 0.156058987 0.001
## 3 BAOR 0.45535259 -0.100215702 0.001
## 4 BARS 0.41614862 0.128649059 0.001
## 5 BCCH -0.25914225 0.186621339 0.001
## 6 BEKI 0.07623645 0.304260844 0.001
## 7 CARW 0.19634245 -0.358991916 0.001
## 8 CEDW 0.50793156 -0.190529183 0.001
## 9 COGR 0.41281036 0.086914877 0.001
## 10 COYE 0.54462857 0.055690572 0.001
## 11 DEJU -0.37274656 0.195628876 0.001
## 12 EAKI 0.38558309 0.275585839 0.001
## 13 EAWP -0.54367957 0.191540154 0.001
## 14 EUST 0.61024663 0.132116554 0.001
## 15 GRCA 0.53815426 -0.159427712 0.001
## 16 HAWO -0.30118219 0.004202757 0.001
## 17 HOFI 0.29177592 0.044619987 0.001
## 18 HOWA -0.44401583 0.047422992 0.001
## 19 HOWR 0.28559364 -0.072630263 0.001
## 20 INBU 0.35169318 -0.145871019 0.001
## 21 KILL 0.25517824 0.354983330 0.001
## 22 MODO 0.51137295 0.027654563 0.001
## 23 NOFL 0.26754080 0.099347369 0.001
## 24 NRWS 0.32290594 0.402369183 0.001
## 25 OROR 0.66614724 -0.071544010 0.001
## 26 PIWO -0.19455428 0.313914955 0.001
## 27 PUFI 0.24916888 0.638570907 0.001
## 28 RBGR 0.31544337 -0.032567502 0.001
## 29 REVI -0.40938637 0.140627611 0.001
## 30 RWBL 0.62715213 0.018320407 0.001
## 31 SCTA -0.55217353 0.233164863 0.001
## 32 SOSP 0.60068008 -0.143602772 0.001
## 33 TRES 0.50591133 0.231005780 0.001
## 34 WAVI 0.63303103 0.006395845 0.001
## 35 WBNU -0.48857239 -0.043518576 0.001
## 36 WIFL 0.49797600 0.196922661 0.001
## 37 WOTH -0.31149124 -0.067476593 0.001
## 38 YEWA 0.72047250 0.049710859 0.001fit_community <- community_fit %>%
scores(., display = "factors") %>%
as.data.frame() %>%
rownames_to_column("community")
rda_communityn_plot_new <- ggplot(plot_df_community, aes(x = RDA1, y = RDA2)) +
coord_fixed() +
geom_point(aes(color = community), size = 3, alpha = 0.8) +
geom_segment(data = fit_community_spp, aes(x = 0, xend = RDA1, y = 0, yend = RDA2),
col = "black") +
geom_text(data = fit_community_spp, aes(label = species))+
geom_text(data=as.data.frame(fit_community),aes(label=community))
rda_communityn_plot_new
ggplot(plot_df_community, aes(x = RDA1, y = RDA2)) +
coord_fixed() +
geom_text(data=as.data.frame(fit_community),aes(label=community))+xlim(-2,4)+ylim(-1,5)
ggplot(plot_df_community, aes(x = RDA1, y = RDA2)) +
coord_fixed() +
geom_text(data=as.data.frame(fit_community_spp),aes(label=species))+xlim(-1,1)+ylim(-1,1)
For these ordinations, RDA1 seems to differentiate birds along an approximate gradient of forest cover - with negative loadings indicating more hardwood-associated bird communities and more positive loadings indicating more mesic/meadow-associated areas. Intriguingly, RDA2 then groups Mesic Meadows out from all other community types - with more positive values indicating more Mesic Meadow-associated birds.
Neotropical migrants: Acadian Fly-Catcher, Baltimore Oriole, Barn Swallow, Common Yellowthroat, Eastern Kingbird, Eastern Wood-Pewee, Hooded Warbler, Indigo Bunting, Orchard Oriole, Rose-Breasted Grosbeak, Red-Eyed Vireo, Scarlet Tanager, Warbling Vireo, Willow Flycatcher, Wood Thrush, Yellow Warbler = 16 species
Resident: Carolina Wren, European Starling, Hairy Woodpecker, House Sparrow, Pileated Woodpecker, Song Sparrow, White-Breasted Nuthatch = 7 species
Short-distance migrant: Brown Thrasher, Cedar Waxwing, Common Grackle, Gray Catbird, House Wren, Killdeer, Mourning Dove, Northern Rough-winged Swallow, Purple Finch, Red-winged Blackbird, Tree Swallow = 11 species
Resident or migratory: Dark-eyed Junco = 1 species
We also can explore linear discriminant analysis:
community<-factor(CMP_bbs_pcap_join_pcaponly_pcapsum$community.x)
birdcomm_bin_drop1<-birdcomm_bin[colSums(birdcomm_bin_min3) > 1]
color_resid_drop1<-color_resid[(colnames(birdcomm_bin) %in% colnames(birdcomm_bin_drop1))]
lda_birdcomm<-lda(birdcomm_bin_drop1, grouping = community)
lda_birdcomm_full<-lda(birdcomm_bin, grouping = community)
ggord(lda_birdcomm, community, xlim=c(-15,30),ylim = c(-15, 15),txt=NULL,arrow=NULL)
ggord(lda_birdcomm, community, xlim=c(-15,30),ylim = c(-15, 15),txt=NULL,arrow=NULL,alpha=0,alpha_el=0.5)
ggord(lda_birdcomm, community, xlim=c(-15,30),ylim = c(-15, 15),alpha=0,alpha_el=0.5,labcol=color_resid_drop1)
ggord(lda_birdcomm_full, community, xlim=c(-15,30),ylim = c(-15, 15),alpha=0,alpha_el=0.5,labcol=color_resid)
ggord(lda_birdcomm_full, community, xlim=c(-10,25),ylim = c(-7.5, 12.5),alpha=0,alpha_el=0.5,var_sub = c("ACFL","HOWA","PIWO","SOSP"))
co-occurrence networks
We can also see which bird species co-occur together.
CMP_bbs_pcap_join_pcaponly_commsumm<-CMP_bbs_pcap_join_pcaponly %>%
group_by(community,species_4_) %>%
summarise(n=n()) %>%
spread(species_4_,n) %>%
replace(is.na(.),0)
CMP_bbs_pcap_join_pcaponly_commsumm<-column_to_rownames(CMP_bbs_pcap_join_pcaponly_commsumm,var = "community")
CMP_bbs_pcap_join_pcaponly_commsumm<-CMP_bbs_pcap_join_pcaponly_commsumm[,-1]
CMP_bbs_pcap_join_pcaponly_commsumm_presabs<-CMP_bbs_pcap_join_pcaponly_commsumm %>% mutate_if(is.numeric, ~1 * (. >= 1))
CMP_bbs_pcap_join_pcaponly_commsumm_presabs_t<-t(CMP_bbs_pcap_join_pcaponly_commsumm_presabs)
CMP_bbs_pcap_join_pcaponly_reservsumm<-CMP_bbs_pcap_join_pcaponly %>%
group_by(reservatio,species_4_) %>%
summarise(n=n()) %>%
spread(species_4_,n) %>%
replace(is.na(.),0)
CMP_bbs_pcap_join_pcaponly_reservsumm<-column_to_rownames(CMP_bbs_pcap_join_pcaponly_reservsumm,var = "reservatio")
CMP_bbs_pcap_join_pcaponly_reservsumm<-CMP_bbs_pcap_join_pcaponly_reservsumm[,-1]
CMP_bbs_pcap_join_pcaponly_reservsumm_presabs<-CMP_bbs_pcap_join_pcaponly_reservsumm %>% mutate_if(is.numeric, ~1 * (. >= 1))
CMP_bbs_pcap_join_pcaponly_reservsumm_presabs_t<-t(CMP_bbs_pcap_join_pcaponly_reservsumm_presabs)
birdcomm_bin_t<-t(birdcomm_bin)
BBSMigrants<-read.csv("G:/NaturalResources/Byer/datasets/BBS/BBS Migants.csv",header=T)
BBSMigrants_sub<-BBSMigrants[BBSMigrants$SpeciesCode %in% colnames(birdcomm_bin),]
BBSMigrants_sub_reorder<-BBSMigrants_sub[match(rownames(birdcomm_bin_t), BBSMigrants_sub$SpeciesCode),]
Residency<-BBSMigrants_sub_reorder$Residency
color_resid<-Residency
color_resid[color_resid %in% "NE"]<-"blue"
color_resid[color_resid %in% "M"]<-"green"
color_resid[color_resid %in% "R"]<-"red"
# presence/absence across plots
co <- print(cooccur(birdcomm_bin_t, spp_names = TRUE))
co[, "sp1_name"] == rownames(birdcomm_bin_t)[co$sp1]
co[, "sp2_name"] == rownames(birdcomm_bin_t)[co$sp2]
nodes <- data.frame(id = 1:nrow(birdcomm_bin_t),
label = rownames(birdcomm_bin_t),
color = color_resid,
shadow = TRUE)
edges <- data.frame(from = co$sp1, to = co$sp2,
color = ifelse(co$p_lt <= 0.05, "#B0B2C1", "#3C3F51"),
dashes = ifelse(co$p_lt <= 0.05, TRUE, FALSE))
cooc<-cooccur(birdcomm_bin_t, spp_names = TRUE)
coccur
# presence/absence across reservations
co_reserv <- print(cooccur(CMP_bbs_pcap_join_pcaponly_reservsumm_presabs_t, spp_names = TRUE))
co_reserv[, "sp1_name"] == rownames(CMP_bbs_pcap_join_pcaponly_reservsumm_presabs_t)[co_reserv$sp1]
co_reserv[, "sp2_name"] == rownames(CMP_bbs_pcap_join_pcaponly_reservsumm_presabs_t)[co_reserv$sp2]
nodes_reserv <- data.frame(id = 1:nrow(CMP_bbs_pcap_join_pcaponly_reservsumm_presabs_t),
label = rownames(CMP_bbs_pcap_join_pcaponly_reservsumm_presabs_t),
color = color_resid,
shadow = TRUE)
edges_reserv <- data.frame(from = co_reserv$sp1, to = co_reserv$sp2,
color = ifelse(co_reserv$p_lt <= 0.05, "#B0B2C1", "#3C3F51"),
dashes = ifelse(co_reserv$p_lt <= 0.05, TRUE, FALSE))
coocc_reserv<-cooccur(CMP_bbs_pcap_join_pcaponly_reservsumm_presabs_t, spp_names = TRUE)
# presence/absence across habitats
co_comm <- print(cooccur(CMP_bbs_pcap_join_pcaponly_commsumm_presabs_t, spp_names = TRUE))
co_comm[, "sp1_name"] == rownames(CMP_bbs_pcap_join_pcaponly_commsumm_presabs_t)[co_comm$sp1]
co_comm[, "sp2_name"] == rownames(CMP_bbs_pcap_join_pcaponly_commsumm_presabs_t)[co_comm$sp2]
nodes_comm <- data.frame(id = 1:nrow(CMP_bbs_pcap_join_pcaponly_commsumm_presabs_t),
label = rownames(CMP_bbs_pcap_join_pcaponly_commsumm_presabs_t),
color = color_resid,
shadow = TRUE)
edges_comm <- data.frame(from = co_comm$sp1, to = co_comm$sp2,
color = ifelse(co_comm$p_lt <= 0.05, "#B0B2C1", "#3C3F51"),
dashes = ifelse(co_comm$p_lt <= 0.05, TRUE, FALSE))
coocc_comm<-cooccur(CMP_bbs_pcap_join_pcaponly_commsumm_presabs_t, spp_names = TRUE)visNetwork(nodes = nodes_comm, edges = edges_comm) %>%
visIgraphLayout(layout = "layout_with_kk")plot(coocc_comm)
plot(coocc_reserv)
Now pRDA - with variance partitioning.
As a final sanity check, it can sometimes be difficult to disentangle geographic effects from environmental effects in these sorts of analyses. For example, what if (hypothetically) cover of seedless vascular plants is negatively associated with geographic coordinates? Thus, we might want to briefly check if the variance explained by our environmental axes is (or is not) confounded with geography. To do so, we will use the varpart() function, which partitions variance between model components:
latlong<-CMP_bbs_pcap_join_pcaponly_pcapsum[,3:4]
vegecomm_latlong<-CMP_bbs_pcap_join_pcaponly_pcapsum[,3:26]
landvegcomm_latlong<-cbind(landcomm,vegecomm,latlong)
pcap_landveggeog<-varpart(birdcomm_bin,~ vibifq+vibifqnotrees+carex+cyper+dicot+
shade+shrub+hydro+svp+ap+fqi+bryo+hydro_rel+hydro_rel_notrees+sens_rel+
sens_rel_notree+tol_rel+tol_rel_notree+inv+small_tree+subcanopy,~developed+activity.area+
closest.developed+sanctioned.trail+bootleg.trail+closest.trail+developed.edge+stream.edge+
total.edges+fragment.area+fragment.perimeter+ldi_3k,~longitude+latitude,data=landvegcomm_latlong)
plot(pcap_landveggeog)
So - most of this variation is associated with these top covariates (X1), and NOT with latitude and longitude alone (X2). Good to see - this means we can be confident that these are actual environmental effects!
community richness metrics
Next, we can explore the relationships between species richness and plot characteristics. We will start with plot-level explorations.
plot-level analyses
plot-level richness and various indicators
Now, we will produce some exploratory models and plots of richness against basic indicators of vegetative quality, including the Floristic Quality Indicator (FQI) and the Vegetative Index of Biotic Integrity (VIBI).
Somewhat unexpectedly, this seems to suggest that higher VIBI scores actually produce slightly less species-rich communities. What could be causing this? Well, to start - perhaps the number of surveys for each plot may be influencing this.
We can generate a summary of the number of surveys for each plot.
So the number of surveys clearly has an effect on the species richness of the site. This is somewhat difficult to account for, but we will attempt to do so using mixed effects models. I then use Akaike’s Information Criterion (AIC for short) to rank these models. The lowest number is the best, and any model within 2 AIC of the top model is potentially worth interpreting.
## Linear mixed model fit by REML ['lmerMod']
## Formula: S_plotdate ~ bootleg.trail + (1 | simp_plot)
## Data: landcomm_bycount
##
## REML criterion at convergence: 4350.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5308 -0.6861 -0.0906 0.6317 3.4835
##
## Random effects:
## Groups Name Variance Std.Dev.
## simp_plot (Intercept) 6.293 2.509
## Residual 13.636 3.693
## Number of obs: 763, groups: simp_plot, 164
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 13.106390 0.300749 43.579
## bootleg.trail -0.000443 0.001394 -0.318
##
## Correlation of Fixed Effects:
## (Intr)
## bootleg.trl -0.594##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt Res.LL
## Mod16 4 4326.31 0.00 0.75 0.75 -2159.13
## Mod2 4 4328.88 2.57 0.21 0.96 -2160.41
## Mod10 4 4332.53 6.22 0.03 0.99 -2162.24
## Mod3 4 4335.67 9.36 0.01 1.00 -2163.81
## Mod9 4 4341.60 15.30 0.00 1.00 -2166.78
## Mod1 4 4343.63 17.32 0.00 1.00 -2167.79
## Mod13 4 4343.66 17.36 0.00 1.00 -2167.81
## Mod14 4 4343.88 17.58 0.00 1.00 -2167.92
## Mod8 4 4344.70 18.39 0.00 1.00 -2168.32
## Mod29 3 4345.27 18.96 0.00 1.00 -2169.62
## Mod17 4 4346.17 19.86 0.00 1.00 -2169.06
## Mod6 4 4346.62 20.31 0.00 1.00 -2169.28
## Mod12 4 4346.68 20.38 0.00 1.00 -2169.32
## Mod26 4 4348.35 22.05 0.00 1.00 -2170.15
## Mod5 4 4349.06 22.75 0.00 1.00 -2170.50
## Mod15 4 4349.16 22.85 0.00 1.00 -2170.55
## Mod4 4 4349.17 22.86 0.00 1.00 -2170.56
## Mod7 4 4350.72 24.42 0.00 1.00 -2171.34
## Mod11 4 4351.45 25.14 0.00 1.00 -2171.70
## Mod18 4 4352.37 26.06 0.00 1.00 -2172.16
## Mod21 4 4358.03 31.72 0.00 1.00 -2174.99
## Mod22 4 4358.07 31.76 0.00 1.00 -2175.01
## Mod20 4 4358.51 32.20 0.00 1.00 -2175.23
## Mod19 4 4360.01 33.71 0.00 1.00 -2175.98
## Mod25 4 4362.87 36.57 0.00 1.00 -2177.41
## Mod27 4 4363.50 37.19 0.00 1.00 -2177.72
## Mod23 4 4364.82 38.52 0.00 1.00 -2178.39
## Mod28 4 4365.53 39.22 0.00 1.00 -2178.74
## Mod24 4 4374.94 48.63 0.00 1.00 -2183.44So richness - at a plot level, and after accounting for plot identity - is determined by how many tolerant species are in the associated plot. And this is much more predictive than any other model.
Here is the plot for this top model:
## `geom_smooth()` using formula = 'y ~ x'
We can also verify that our mixed effects model fits better compared to the fixed effects model.
Rich_tol_rel_fixed<-lm(S_plotdate~tol_rel,vegecomm_bycount)
anova(Rich_tol_rel, Rich_tol_rel_fixed) ## refitting model(s) with ML (instead of REML)## Data: vegecomm_bycount
## Models:
## Rich_tol_rel_fixed: S_plotdate ~ tol_rel
## Rich_tol_rel: S_plotdate ~ tol_rel + (1 | simp_plot)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## Rich_tol_rel_fixed 3 4391.9 4405.8 -2193.0 4385.9
## Rich_tol_rel 4 4326.6 4345.2 -2159.3 4318.6 67.289 1 2.345e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1This confirms that we chose correctly - the mixed effects model fits better than the fixed effects-only model.
We can also inspect how this observation-level mixed effect model compares to a plot-level, fixed effect model (in other words - a view that just looks at richness at a plot level - aggregating over observations).
S_tol_rel<-lm(S~tol_rel,vegecomm_richness)
confint(S_tol_rel)## 2.5 % 97.5 %
## (Intercept) 22.88094 25.853189
## tol_rel 1.96969 9.053297confint(Rich_tol_rel)## Computing profile confidence intervals ...## 2.5 % 97.5 %
## .sig01 1.845565 2.700243
## .sigma 3.493475 3.916023
## (Intercept) 11.100919 12.494418
## tol_rel 2.204396 5.564987We can see that, while the covariate estimate for the plot-level model is higher, it is also estimated with much more variance (its confidence intervals are wider). This is not surprising - our functional sample size for this dataset is 164 plots. In contrast, while the observation-level mixed-effects model has a lower covariate estimate, we have tighter confidence intervals around this estimate - likely by virtue of a larger number of samples and a more effectively parameterized model.
reservation-level analyses
Now we are going to re-summarize this dataset at a reservation-level. While we will have limited power for many comparisons of interest, this could provide some indication of whether covariate relationships vary across scales.
This doesn’t quite get at what we want - this displays species accumulation across reservations after accounting for survey effort per reservation. Let’s try exploring reservation-by-reservation rarefaction with the previous observation-level dataset.
With the plot-level data, we could easily account for variation in survey effort across plots with a mixed model design, but this case is a bit more complicated. We are interested in reservation-level species richness, but know that reservations were visited different numbers of times - based on the number of plots and the number of surveys/plot. Thus, the easiest thing to start with is simply using species richness estimates interpolated to the least number of surveys per reservation (3, in our case).
Now we can proceed with exploring richness patterns across reservations.
##
## Model selection based on AICc:
##
## K AICc Delta_AICc AICcWt Cum.Wt LL
## Mod18 3 80.78 0.00 0.39 0.39 -36.39
## Mod8 3 82.73 1.95 0.15 0.54 -37.36
## Mod16 3 82.82 2.04 0.14 0.68 -37.41
## Mod17 3 83.59 2.80 0.10 0.77 -37.79
## Mod7 3 84.81 4.03 0.05 0.82 -38.40
## Mod12 3 84.97 4.19 0.05 0.87 -38.49
## Mod11 3 87.00 6.22 0.02 0.89 -39.50
## Mod19 2 87.17 6.38 0.02 0.91 -41.12
## Mod1 3 87.38 6.60 0.01 0.92 -39.69
## Mod6 3 87.62 6.84 0.01 0.93 -39.81
## Mod10 3 87.69 6.91 0.01 0.95 -39.85
## Mod14 3 87.94 7.16 0.01 0.96 -39.97
## Mod2 3 88.29 7.50 0.01 0.97 -40.14
## Mod4 3 88.85 8.06 0.01 0.97 -40.42
## Mod15 3 88.85 8.07 0.01 0.98 -40.43
## Mod13 3 89.25 8.47 0.01 0.99 -40.63
## Mod5 3 89.30 8.52 0.01 0.99 -40.65
## Mod3 3 89.50 8.71 0.00 1.00 -40.75
## Mod9 3 89.83 9.04 0.00 1.00 -40.91
although our power is a bit limited above, it does appear that the Vegetative Index of Biotic Integrity - at a reservation-level - does appear to predict bird species richness (as long as we do not include woody plants in that list!)
we can also plot out richness for each site:
In confirmation of our prior ordination plots, it now appears that Rocky River is one of the most species rich areas - followed by Hinckley, South Chagrin, Bedford, and Brecksville.
multi-species occupancy modelling
MSOM - running different models
setting some basic parameters for the runs
N <- dim(bbs_occobject$y)[1]
ms.inits <- list(alpha.comm = 0,
beta.comm = 0,
beta = 0,
alpha = 0,
tau.sq.beta = 1,
tau.sq.alpha = 1,
z = apply(bbs_occobject$y, c(1, 2), max, na.rm = TRUE))
ms.priors <- list(beta.comm.normal = list(mean = 0, var = 2.72),
alpha.comm.normal = list(mean = 0, var = 2.72),
tau.sq.beta.ig = list(a = 0.1, b = 0.1),
tau.sq.alpha.ig = list(a = 0.1, b = 0.1))
n.samples<-30000
n.omp.threads<-1
n.report<-6000
n.burn<-10000
n.thin<-50
n.chains<-3setting up our model statements
occ_null <- ~1
occ_vibifq <- ~ scale(vibifq)
occ_vibifqnotrees <- ~ scale(vibifqnotrees)
occ_fqi <- ~ scale(fqi)
occ_rdaveg <- ~ scale(tol_rel) + scale(svp) + scale(dicot) + scale(cyper) + scale(small_tree)
occ_fullrda<- ~ scale(tol_rel) + scale(svp) + scale(dicot) + scale(cyper) +scale(ldi_3k)+scale(developed)
occ_structure <- ~ scale(canopy)+scale(subcanopy)+scale(small_tree)
occ_tol <- ~ scale(tol_rel)
occ_allveg <- ~ scale(tol_rel) + scale(svp) + scale(dicot) + scale(cyper) + scale(small_tree)+scale(shrub)+scale(hydro)+scale(subcanopy)+scale(canopy)
occ_allland <- ~ scale(developed)+scale(activity.area)+scale(sanctioned.trail)+ scale(bootleg.trail)+scale(stream.edge)+scale(ldi_3k)
det_null <- ~1
det_null_randobs<- ~ 1 + I(1|observer)
det_day <- ~ scale(day)
det_day_randobs <- ~ scale(day) + I(1|observer)
det_daytod<- ~ scale(day)+scale(tod)
det_daytod_randobs <- ~ scale(day)+scale(tod)+ I(1|observer)
det_daytodquad <- ~ scale(day)+scale(tod) +I(scale(day)^2)
det_daytodquad_randobs <- ~ scale(day)+scale(tod) +I(scale(day)^2) + I(1|observer)
det_season <- ~ scale(seasons)
det_season_randobs <- ~ scale(seasons) + I(1|observer)
det_dayyearseason <- ~ scale(dayyear) + scale(seasons)
det_dayyearseason_randobs <- ~ scale(dayyear) + scale(seasons) + I(1|observer)
det_dayyeartodseason <- ~ scale(dayyear)+scale(tod)+scale(seasons)
det_dayyeartodseason_randobs <- ~ scale(dayyear)+scale(tod)+scale(seasons) + I(1|observer)
det_dayyeartodseason <- ~ scale(dayyear)+scale(tod)+scale(seasons)
det_dayyeartodseason_randobs <- ~ scale(dayyear)+scale(tod)+scale(seasons)+ I(1|observer)
det_dayyeartodseasonroad <- ~ scale(dayyear)+scale(tod)+scale(seasons)+scale(dist_road)
det_dayyeartodseasonroad_randobs <- ~ scale(dayyear)+scale(tod)+scale(seasons)+scale(dist_road)+I(1|observer)basic null model
out.ms.null <- msPGOcc(occ.formula = occ_null,
det.formula = det_null,
data = bbs_occobject,
inits = ms.inits,
n.samples = n.samples,
priors = ms.priors,
n.omp.threads = n.omp.threads,
verbose = TRUE,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)
summary(out.ms.null)
ppc.out.ms.null <- ppcOcc(out.ms.null, 'freeman-tukey', group = 1)
summary(ppc.out.ms.null)basic null model - with random effect
out.ms.nullobs <- msPGOcc(occ.formula = occ_null,
det.formula = det_null_randobs,
data = bbs_occobject,
inits = ms.inits,
n.samples = n.samples,
priors = ms.priors,
n.omp.threads = n.omp.threads,
verbose = TRUE,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)
summary(out.ms.nullobs)
ppc.out.ms.nullobs <- ppcOcc(out.ms.nullobs, 'freeman-tukey', group = 1)
summary(ppc.out.ms.nullobs)building up complexity of detection model - just day for now
out.ms.det.day <- msPGOcc(occ.formula = occ_null,
det.formula = det_day,
data = bbs_occobject,
inits = ms.inits,
n.samples = n.samples,
priors = ms.priors,
n.omp.threads = n.omp.threads,
verbose = TRUE,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)
summary(out.ms.det.day)
ppc.out.ms.det.day <- ppcOcc(out.ms.det.day, 'freeman-tukey', group = 1)
summary(ppc.out.ms.det.day)day and random effect
out.ms.detobs.day <- msPGOcc(occ.formula = occ_null,
det.formula = det_day_randobs,
data = bbs_occobject,
inits = ms.inits,
n.samples = n.samples,
priors = ms.priors,
n.omp.threads = n.omp.threads,
verbose = TRUE,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)
summary(out.ms.detobs.day)
ppc.out.ms.detobs.day <- ppcOcc(out.ms.detobs.day, 'freeman-tukey', group = 1)
summary(ppc.out.ms.detobs.day)day and time-of-day
out.ms.det.daytod <- msPGOcc(occ.formula = occ_null,
det.formula = det_daytod,
data = bbs_occobject,
inits = ms.inits,
n.samples = n.samples,
priors = ms.priors,
n.omp.threads = n.omp.threads,
verbose = TRUE,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)
summary(out.ms.det.daytod)
ppc.out.ms.det.daytod <- ppcOcc(out.ms.det.daytod, 'freeman-tukey', group = 1)
summary(ppc.out.ms.det.daytod)day and time-of-day - random observer
out.ms.detobs.daytod <- msPGOcc(occ.formula = occ_null,
det.formula = det_daytod_randobs,
data = bbs_occobject,
inits = ms.inits,
n.samples = n.samples,
priors = ms.priors,
n.omp.threads = n.omp.threads,
verbose = TRUE,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)
summary(out.ms.detobs.daytod)
ppc.out.ms.detobs.daytod <- ppcOcc(out.ms.detobs.daytod, 'freeman-tukey', group = 1)
summary(ppc.out.ms.detobs.daytod)day, tod, quadratic day
out.ms.det.daytodquad<- msPGOcc(occ.formula = occ_null,
det.formula = det_daytodquad,
data = bbs_occobject,
inits = ms.inits,
n.samples = n.samples,
priors = ms.priors,
n.omp.threads = n.omp.threads,
verbose = TRUE,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)
summary(out.ms.det.daytodquad)
ppc.out.ms.det.daytodquad <- ppcOcc(out.ms.det.daytodquad, 'freeman-tukey', group = 1)
summary(ppc.out.ms.det.daytodquad)day, tod, quadratic day - random observer
out.ms.detobs.daytodquad<- msPGOcc(occ.formula = occ_null,
det.formula = det_daytodquad_randobs,
data = bbs_occobject,
inits = ms.inits,
n.samples = n.samples,
priors = ms.priors,
n.omp.threads = n.omp.threads,
verbose = TRUE,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)
summary(out.ms.detobs.daytodquad)
ppc.out.ms.detobs.daytodquad <- ppcOcc(out.ms.detobs.daytodquad, 'freeman-tukey', group = 1)
summary(ppc.out.ms.detobs.daytodquad)season alone
out.ms.det.season<- msPGOcc(occ.formula = occ_null,
det.formula = det_season,
data = bbs_occobject,
inits = ms.inits,
n.samples = n.samples,
priors = ms.priors,
n.omp.threads = n.omp.threads,
verbose = TRUE,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)
summary(out.ms.det.season)
ppc.out.ms.det.season <- ppcOcc(out.ms.det.season, 'freeman-tukey', group = 1)
summary(ppc.out.ms.det.season)season alone - random observer
out.ms.detobs.season<- msPGOcc(occ.formula = occ_null,
det.formula = det_season_randobs,
data = bbs_occobject,
inits = ms.inits,
n.samples = n.samples,
priors = ms.priors,
n.omp.threads = n.omp.threads,
verbose = TRUE,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)
summary(out.ms.detobs.season)
ppc.out.ms.detobs.season <- ppcOcc(out.ms.detobs.season, 'freeman-tukey', group = 1)
summary(ppc.out.ms.detobs.season)day of year and season
out.ms.det.dayyearseason<- msPGOcc(occ.formula = occ_null,
det.formula = det_dayyearseason,
data = bbs_occobject,
inits = ms.inits,
n.samples = n.samples,
priors = ms.priors,
n.omp.threads = n.omp.threads,
verbose = TRUE,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)
summary(out.ms.det.dayyearseason)
ppc.out.ms.det.dayyearseason <- ppcOcc(out.ms.det.dayyearseason, 'freeman-tukey', group = 1)
summary(ppc.out.ms.det.dayyearseason)day of year and season - random observer
out.ms.detobs.dayyearseason<- msPGOcc(occ.formula = occ_null,
det.formula = det_dayyearseason_randobs,
data = bbs_occobject,
inits = ms.inits,
n.samples = n.samples,
priors = ms.priors,
n.omp.threads = n.omp.threads,
verbose = TRUE,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)
summary(out.ms.detobs.dayyearseason)
ppc.out.ms.detobs.dayyearseason <- ppcOcc(out.ms.detobs.dayyearseason, 'freeman-tukey', group = 1)
summary(ppc.out.ms.detobs.dayyearseason)day in year, tod, season
out.ms.det.dayyeartodseason<- msPGOcc(occ.formula = occ_null,
det.formula = det_dayyeartodseason,
data = bbs_occobject,
inits = ms.inits,
n.samples = n.samples,
priors = ms.priors,
n.omp.threads = n.omp.threads,
verbose = TRUE,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)
summary(out.ms.det.dayyeartodseason)
ppc.out.ms.det.dayyeartodseason <- ppcOcc(out.ms.det.dayyeartodseason, 'freeman-tukey', group = 1)
summary(ppc.out.ms.det.dayyeartodseason)day in year, tod, season - random observer
out.ms.detobs.dayyeartodseason<- msPGOcc(occ.formula = occ_null,
det.formula = det_dayyeartodseason_randobs,
data = bbs_occobject,
inits = ms.inits,
n.samples = n.samples,
priors = ms.priors,
n.omp.threads = n.omp.threads,
verbose = TRUE,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)
summary(out.ms.detobs.dayyeartodseason)
ppc.out.ms.detobs.dayyeartodseason <- ppcOcc(out.ms.detobs.dayyeartodseason, 'freeman-tukey', group = 1)
summary(ppc.out.ms.detobs.dayyeartodseason)day in year, tod, season
out.ms.det.dayyeartodseasonroad<- msPGOcc(occ.formula = occ_null,
det.formula = det_dayyeartodseasonroad,
data = bbs_occobject,
inits = ms.inits,
n.samples = n.samples,
priors = ms.priors,
n.omp.threads = n.omp.threads,
verbose = TRUE,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)
summary(out.ms.det.dayyeartodseasonroad)
ppc.out.ms.det.dayyeartodseasonroad <- ppcOcc(out.ms.det.dayyeartodseasonroad, 'freeman-tukey', group = 1)
summary(ppc.out.ms.det.dayyeartodseasonroad)day in year, tod, season - random observer
out.ms.detobs.dayyeartodseasonroad<- msPGOcc(occ.formula = occ_null,
det.formula = det_dayyeartodseasonroad_randobs,
data = bbs_occobject,
inits = ms.inits,
n.samples = n.samples,
priors = ms.priors,
n.omp.threads = n.omp.threads,
verbose = TRUE,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)
summary(out.ms.detobs.dayyeartodseasonroad)
ppc.out.ms.detobs.dayyeartodseasonroad <- ppcOcc(out.ms.detobs.dayyeartodseasonroad, 'freeman-tukey', group = 1)
summary(ppc.out.ms.detobs.dayyeartodseasonroad)selecting fixed effects detection model
waicOcc(out.ms.null)
waicOcc(out.ms.det.day)
waicOcc(out.ms.det.daytod)
waicOcc(out.ms.det.daytodquad)
waicOcc(out.ms.det.season)
waicOcc(out.ms.det.dayyearseason)
waicOcc(out.ms.det.dayyeartodseason)
waicOcc(out.ms.det.dayyeartodseasonroad)selecting random effect detection model
waicOcc(out.ms.nullobs)
waicOcc(out.ms.detobs.day)
waicOcc(out.ms.detobs.daytod)
waicOcc(out.ms.detobs.daytodquad)
waicOcc(out.ms.detobs.season)
waicOcc(out.ms.detobs.dayyearseason)
waicOcc(out.ms.detobs.dayyeartodseason)
waicOcc(out.ms.detobs.dayyeartodseasonroad)The best-supported model is day of year, time of day, season, and distance to road. Thus, we will proceed with that one for now.
Now we will run several different occupancy submodels, once again evaluating fit for progressively more complicated models.
occupancy model with vibi-fq
det_best<-det_dayyeartodseasonroad_randobs
out.ms.vibifq <- msPGOcc(occ.formula = occ_vibifq,
det.formula = det_best,
data = bbs_occobject,
inits = ms.inits,
n.samples = n.samples,
priors = ms.priors,
n.omp.threads = n.omp.threads,
verbose = TRUE,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)
summary(out.ms.vibifq)
ppc.out.ms.vibifq <- ppcOcc(out.ms.vibifq, 'freeman-tukey', group = 1)
summary(ppc.out.ms.vibifq)
waicOcc(out.ms.vibifq)occupancy model with vibi-fq-notrees
out.ms.vibifqnotrees <- msPGOcc(occ.formula = occ_vibifqnotrees,
det.formula = det_best,
data = bbs_occobject,
inits = ms.inits,
n.samples = n.samples,
priors = ms.priors,
n.omp.threads = n.omp.threads,
verbose = TRUE,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)
summary(out.ms.vibifqnotrees)
ppc.out.ms.vibifqnotrees <- ppcOcc(out.ms.vibifqnotrees, 'freeman-tukey', group = 1)
summary(ppc.out.ms.vibifqnotrees)
waicOcc(out.ms.vibifqnotrees)just FQAI
out.ms.fqi <- msPGOcc(occ.formula = occ_fqi,
det.formula = det_best,
data = bbs_occobject,
inits = ms.inits,
n.samples = n.samples,
priors = ms.priors,
n.omp.threads = n.omp.threads,
verbose = TRUE,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)
summary(out.ms.fqi)
ppc.out.ms.fqi <- ppcOcc(out.ms.fqi, 'freeman-tukey', group = 1)
summary(ppc.out.ms.fqi)
waicOcc(out.ms.fqi)occupancy model with all rda-supported vegetation covariates
out.ms.rdaveg <- msPGOcc(occ.formula = occ_rdaveg,
det.formula = det_best,
data = bbs_occobject,
inits = ms.inits,
n.samples = n.samples,
priors = ms.priors,
n.omp.threads = n.omp.threads,
verbose = TRUE,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)
summary(out.ms.rdaveg)
ppc.out.ms.rdaveg <- ppcOcc(out.ms.rdaveg, 'freeman-tukey', group = 1)
summary(ppc.out.ms.rdaveg)
waicOcc(out.ms.rdaveg)occupancy model with all rda-supported covariates
out.ms.fulllandveg <- msPGOcc(occ.formula = occ_fullrda,
det.formula = det_best,
data = bbs_occobject,
inits = ms.inits,
n.samples = n.samples,
priors = ms.priors,
n.omp.threads = n.omp.threads,
verbose = TRUE,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)
summary(out.ms.fulllandveg)
ppc.out.ms.fulllandveg <- ppcOcc(out.ms.fulllandveg, 'freeman-tukey', group = 1)
summary(ppc.out.ms.fulllandveg)
waicOcc(out.ms.fulllandveg)
write_rds(out.ms.fulllandveg,"G:/NaturalResources/Byer/datasets/BBS/out.ms.fulllandveg.rds")occupancy model with just tolerance/sensitivity
out.ms.tolerance<- msPGOcc(occ.formula = occ_tol,
det.formula = det_best,
data = bbs_occobject,
inits = ms.inits,
n.samples = n.samples,
priors = ms.priors,
n.omp.threads = n.omp.threads,
verbose = TRUE,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)
summary(out.ms.tolerance)
ppc.out.ms.tolerance <- ppcOcc(out.ms.tolerance, 'freeman-tukey', group = 1)
summary(ppc.out.ms.tolerance)
waicOcc(out.ms.tolerance)all veg covariates model
out.ms.allveg <- msPGOcc(occ.formula = occ_allveg,
det.formula = det_best,
data = bbs_occobject,
inits = ms.inits,
n.samples = n.samples,
priors = ms.priors,
n.omp.threads = n.omp.threads,
verbose = TRUE,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)
summary(out.ms.allveg)
ppc.out.ms.allveg <- ppcOcc(out.ms.allveg, 'freeman-tukey', group = 1)
summary(ppc.out.ms.allveg)
waicOcc(out.ms.allveg)all landscape covariates
out.ms.land <- msPGOcc(occ.formula = occ_allland,
det.formula = det_best,
data = bbs_occobject,
inits = ms.inits,
n.samples = n.samples,
priors = ms.priors,
n.omp.threads = n.omp.threads,
verbose = TRUE,
n.report = n.report,
n.burn = n.burn,
n.thin = n.thin,
n.chains = n.chains)
summary(out.ms.land)
ppc.out.ms.land <- ppcOcc(out.ms.land, 'freeman-tukey', group = 1)
summary(ppc.out.ms.land)
waicOcc(out.ms.land)model selection
we can also print out a table of Widely Applicable Information Criterion (wAIC) values. The lower this is, the better the model fit.
AICc_table<-tibble(.rows = 9)
AICc_table$model<-c("Null", "p ~ day + tod + season + dist_roads, psi ~ tol_rel + svp + dicot + cyper + small_tree + shrub + hydro + subcanopy + canopy", "p ~ day + tod + season + dist_roads, psi ~ tol_rel + svp + dicot + cyper + ldi + developed","p ~ day + tod + season + dist_roads, psi ~ developed + activityarea + sanctionedtrail + bootlegtrail + streamedge + ldi","p ~ day + tod + season + dist_roads,psi~tol_rel + svp + dicot + cyper + small_tree","p ~ day + tod + season + dist_roads,psi~canopy + subcanopy + small_tree","p ~ day + tod + season + dist_roads,psi~tol_rel","p ~ day + tod + season + dist_roads,psi~vibi","p ~ day + tod + season + dist_roads,psi~vibi_notrees")
AICc_table$wAIC<-c(waicOcc(out.ms.nullobs)[[3]],waicOcc(out.ms.allveg)[[3]],waicOcc(out.ms.fulllandveg)[[3]],waicOcc(out.ms.land)[[3]],waicOcc(out.ms.rdaveg)[[3]],
waicOcc(out.ms.structure)[[3]],waicOcc(out.ms.tolerance)[[3]],waicOcc(out.ms.vibifq)[[3]],waicOcc(out.ms.vibifqnotrees)[[3]])| model | wAIC | |
|---|---|---|
| 3 | p ~ day + tod + season + dist_roads, psi ~ tol_rel + svp + dicot + cyper + ldi + developed | 35681.75 |
| 2 | p ~ day + tod + season + dist_roads, psi ~ tol_rel + svp + dicot + cyper + small_tree + shrub + hydro + subcanopy + canopy | 35854.36 |
| 5 | p ~ day + tod + season + dist_roads,psi~tol_rel + svp + dicot + cyper + small_tree | 35874.89 |
| 8 | p ~ day + tod + season + dist_roads,psi~vibi | 35882.12 |
| 7 | p ~ day + tod + season + dist_roads,psi~tol_rel | 36040.55 |
| 9 | p ~ day + tod + season + dist_roads,psi~vibi_notrees | 36066.59 |
| 4 | p ~ day + tod + season + dist_roads, psi ~ developed + activityarea + sanctionedtrail + bootlegtrail + streamedge + ldi | 36172.77 |
| 6 | p ~ day + tod + season + dist_roads,psi~canopy + subcanopy + small_tree | 36402.48 |
| 1 | Null | 36879.23 |
So, the top model is the model with a mix of land cover and vegetation covariates!
Now, we plot based on this top model - after running it for a bit longer.
out.ms.fulllandveg_2 <- msPGOcc(occ.formula = occ_fullrda,
det.formula = det_best,
data = bbs_occobject,
inits = ms.inits,
n.samples = 200000,
priors = ms.priors,
n.omp.threads = 4,
verbose = TRUE,
n.report = 190000,
n.burn = 10000,
n.thin = 100,
n.chains = n.chains)
summary(out.ms.fulllandveg_2)
ppc.out.ms.fulllandveg_2 <- ppcOcc(out.ms.fulllandveg_2, 'freeman-tukey', group = 1)
summary(ppc.out.ms.fulllandveg_2)
waicOcc(out.ms.fulllandveg_2)
write_rds(out.ms.fulllandveg_2,"G:/NaturalResources/Byer/datasets/BBS/out.ms.fulllandveg_2.rds")This top model still has some diagnostics of adequate, but marginal fit, so we will run this without the random effect - and for fewer surveys.
sitechecks_road_detcov_three<-sitechecks_road_detcov[,1:3]
bbs_occobject_three$det.covs$dist_road<-sitechecks_road_detcov_three
out.ms.fulllandveg_three <- msPGOcc(occ.formula = occ_fullrda,
det.formula = det_dayyeartodseasonroad,
data = bbs_occobject_three,
inits = ms.inits,
n.samples = 20000,
priors = ms.priors,
n.omp.threads = 1,
verbose = TRUE,
n.report = 19000,
n.burn = 1000,
n.thin = 50,
n.chains = n.chains)
summary(out.ms.fulllandveg_three)
ppc.out.ms.fulllandveg_three <- ppcOcc(out.ms.fulllandveg_three, 'freeman-tukey', group = 1)
summary(ppc.out.ms.fulllandveg_three)
waicOcc(out.ms.fulllandveg_three)
write_rds(out.ms.fulllandveg_three,"G:/NaturalResources/Byer/datasets/BBS/out.ms.fulllandveg_three.rds")We can also run a model with the “no tree” analog of tol_rel.
occ_fullrda_notree<-~scale(tol_rel_notree) + scale(svp) + scale(dicot) + scale(cyper) +
scale(ldi_3k) + scale(developed)
out.ms.fulllandveg_three_notree <- msPGOcc(occ.formula = occ_fullrda_notree,
det.formula = det_dayyeartodseasonroad,
data = bbs_occobject_three,
inits = ms.inits,
n.samples = 20000,
priors = ms.priors,
n.omp.threads = 1,
verbose = TRUE,
n.report = 1000,
n.burn = 1000,
n.thin = 50,
n.chains = n.chains)
summary(out.ms.fulllandveg_three_notree)
ppc.out.ms.fulllandveg_three_notree <- ppcOcc(out.ms.fulllandveg_three_notree, 'freeman-tukey', group = 1)
summary(ppc.out.ms.fulllandveg_three_notree)
waicOcc(out.ms.fulllandveg_three_notree)
write_rds(out.ms.fulllandveg_three_notree,"G:/NaturalResources/Byer/datasets/BBS/out.ms.fulllandveg_three_notree.rds")This did not effect overall community-level estimates, and we are likely safe to proceed with the original model.
plot single covariate effects
prediction overall all PCAP
pcap_full<-read.csv("G:/NaturalResources/Byer/datasets/BBS/pcap_2nd_sample_database.csv",header=T)
names(pcap_full)[names(pcap_full) == 'Plot.ID'] <- 'simp_plot'
pcap_full_land<-left_join(pcap_full,pcap_1styear_landscape,by="simp_plot")
tol_rel = ((pcap_full_land$tolerant_rel_cov_metric_value - mean(bbs_occobject$occ.covs[,17]))/sd(bbs_occobject$occ.covs[,17]))
svp = ((pcap_full_land$svp_metric_value - mean(bbs_occobject$occ.covs[,9]))/sd(bbs_occobject$occ.covs[,9]))
dicot = ((pcap_full_land$dicot_metric_value - mean(bbs_occobject$occ.covs[, 5])) / sd(bbs_occobject$occ.covs[, 5]))
cyper = ((pcap_full_land$cyperaceae_metric_value - mean(bbs_occobject$occ.covs[, 4])) / sd(bbs_occobject$occ.covs[, 4]))
ldi_3k = ((as.numeric(pcap_full_land$ldi_3k) - mean(bbs_occobject$occ.covs[, 34])) / sd(bbs_occobject$occ.covs[, 34]))
developed = ((as.numeric(pcap_full_land$developed) - mean(bbs_occobject$occ.covs[, 23])) / sd(bbs_occobject$occ.covs[, 23]))
predict_data<-data.frame(int = 1, tol_rel = tol_rel, svp = svp, dicot = dicot, cyper = cyper, ldi_3k = ldi_3k, developed = developed)
fullpcappred<-predict(out.ms.fulllandveg_2,as.matrix(predict_data))
RMSE <- function(observed, predicted) {
sqrt(mean((predicted - observed)^2, na.rm=TRUE))
}
f1 <- function(x, test, train) {
nmx <- x[1]
idp <- x[2]
if (nmx < 1) return(Inf)
if (idp < .001) return(Inf)
m <- gstat(formula=mean.psi~1, locations=train, nmax=nmx, set=list(idp=idp))
p <- predict(m, newdata=test, debug.level=0)$var1.pred
RMSE(test$mean.psi, p)
}
set.seed(12345)
respoly<-readOGR("G:/NaturalResources/Byer/DataAnalysis/bbs_cmp_analyses/bbsapp/CMPboundaries.shp")
fullpcaprasters<-list()
for(i in 1:length(sppnames)){
plot.dat.ms <- data.frame(x = pcap_full_land$xcoord,
y = pcap_full_land$ycoord,
mean.psi = apply(fullpcappred$psi.0.samples[,i,], 2, mean),
sd.psi = apply(fullpcappred$psi.0.samples[,i,], 2, sd),
stringsAsFactors = FALSE)
dsp <- SpatialPoints(plot.dat.ms[,1:2])
dsp <- SpatialPointsDataFrame(dsp, plot.dat.ms)
crs(dsp)<-"+proj=lcc +lat_0=39.6666666666667 +lon_0=-82.5 +lat_1=41.7 +lat_2=40.4333333333333 +x_0=600000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs"
j <- sample(nrow(dsp), 0.2 * nrow(dsp))
tst <- dsp[j,]
trn <- dsp[-j,]
opt <- optim(c(8, .5), f1, test=tst, train=trn)
# define groups for mapping
cuts <- c(0,0.25,0.75,1.0)
# set up a palette of interpolated colors
blues <- colorRampPalette(c('yellow', 'orange', 'blue', 'dark blue'))
spplot(dsp, 'mean.psi', cuts=cuts, col.regions=blues(5), pch=20, cex=2)
v <- voronoi(dsp)
ch <- chull(plot.dat.ms[,1:2])
coords <- plot.dat.ms[c(ch, ch[1]), 1:2] # closed polygon
plot(plot.dat.ms[,1:2], pch=19)
lines(coords, col="red")
#sp_poly <- SpatialPolygons(list(Polygons(list(Polygon(coords)), ID=1)))
# set coordinate reference system with SpatialPolygons(..., proj4string=CRS(...))
# e.g. CRS("+proj=longlat +datum=WGS84")
#sp_poly_df <- SpatialPolygonsDataFrame(sp_poly, data=data.frame(ID=1))
#writeOGR(sp_poly_df, "chull", layer="chull", driver="ESRI Shapefile")
vca <- intersect(v, respoly)
spplot(vca, 'mean.psi', col.regions=rev(get_col_regions()))
r <- raster(respoly, res=500)
crs(r)<-"+proj=lcc +lat_0=39.6666666666667 +lon_0=-82.5 +lat_1=41.7 +lat_2=40.4333333333333 +x_0=600000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs"
crs(vca)<-"+proj=lcc +lat_0=39.6666666666667 +lon_0=-82.5 +lat_1=41.7 +lat_2=40.4333333333333 +x_0=600000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs"
vr <- rasterize(vca, r, 'mean.psi')
crs(vr)<-"+proj=lcc +lat_0=39.6666666666667 +lon_0=-82.5 +lat_1=41.7 +lat_2=40.4333333333333 +x_0=600000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs"
plot(vr)
points(dsp)
writeRaster(vr,filename = paste0("G:/NaturalResources/Byer/datasets/BBS/rasterpredictions/",sppnames[i],"_occmeanpsi_fulllandveg_2_fullpcap.asc"),overwrite=TRUE)
fitmax <- gstat::gstat(formula = mean.psi ~ 1, data = dsp, nmax = opt$par[1], set = list(idp = opt$par[2]))
maxint <- raster::interpolate(r, model=fitmax)
crs(maxint)<-"+proj=lcc +lat_0=39.6666666666667 +lon_0=-82.5 +lat_1=41.7 +lat_2=40.4333333333333 +x_0=600000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs"
plot(maxint, col=rev(heat.colors(255)))
plot(vr, add=TRUE)
r2 <- crop(maxint, extent(respoly))
r3 <- mask(r2, respoly)
fullpcaprasters[[i]]<-r3
writeRaster(r3,filename = paste0("G:/NaturalResources/Byer/datasets/BBS/rasterpredictions/",sppnames[i],"_occmeanpsi_fulllandveg_2_fullpcap_nni_opt_clipped.asc"),overwrite=TRUE)
}
saveRDS(fullpcaprasters,"G:/NaturalResources/Byer/datasets/BBS/fullpcaprasters_occ.Rdata")
#first import all files in a single folder as a list
#rastlist <- list.files(path = "G:/NaturalResources/Byer/datasets/BBS/rasterpredictions/nni_opt/clipped", pattern='.asc$',
#all.files=TRUE, full.names=TRUE)
#import all raster files in folder using lapply
#allrasters <- lapply(rastlist, raster)
#to check the index numbers of all imported raster list elements
#allrasters
#call single raster element
#allrasters[[1]]
#to run a function on an individual raster e.g., plot
#plot(allrasters[[1]])
#saveRDS(allrasters,"G:/NaturalResources/Byer/datasets/BBS/rasterpredictions/nni_opt/clipped/nnirasters_fullpcap.RData")plotting using an RShiny app
While the below RShiny app should be fairly easy to use within the current webpage, you may also access this via: https://testudinidude.shinyapps.io/bbsapp/.
covariate plotting
Finally, we can also run some plots depicting covariate effects for each species. These should make it clear if and when species appear to be particularly related to certain vegetative characteristics - and thus may serve as “indicators” for these characteristics.
just plotting a few focal species
